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ABSTRACT 

A fast algorithm is developed for the directional correlation of scalar band-limited signals and band-limited 
steerable filters on the sphere. The asymptotic complexity associated to it through simple quadrature is of 
order 0{L 5 ), where 2L stands for the square-root of the number of sampling points on the sphere, also setting 
a band limit L for the signals and filters considered. The filter steerability allows to compute the directional 
correlation uniquely in terms of direct and inverse scalar spherical harmonics transforms, which drive the over- 
all asymptotic complexity. The separation of variables technique for the scalar spherical harmonics transform 
produces an 0{L 3 ) algorithm independently of the pixelization. On equi-angular pixelizations, a sampling the- 
orem introduced by Driscoll and Healy implies the exactness of the algorithm. The equi-angular and HEALPix 
implementations are compared in terms of memory requirements, computation times, and numerical stability. 
The computation times for the scalar transform, and hence for the directional correlation, of maps of several 
megapixels on the sphere (L ~ 10 3 ) are reduced from years to tens of seconds in both implementations on a sin- 
gle standard computer. These generic results for the scale-space signal processing on the sphere are specifically 
developed in the perspective of the wavelet analysis of the cosmic microwave background (CMB) temperature 
(T) and polarization (E and B) maps of the WMAP and Planck experiments. As an illustration, we consider the 
computation of the wavelet coefficients of a simulated temperature map of several megapixels with the second 
Gaussian derivative wavelet. 

Subject headings: cosmology: cosmic microwave background — methods: data analysis — methods: numeri- 
cal 
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1. INTRODUCTION 

The cosmic microwave background (CMB) temperature 
and polarization anisotropy distribution represents today an 
exceptional cosmological laboratory for the study of the ge- 
ometrical structure, energy content, and evolution of the uni- 
verse. The last fifteen years of detection and analysis, cul- 
minating with the all-sky maps of the Wilkinson Microwave 
Anisotropy Probe (WMAP) satellite mission, has led to a 
rather precise determination of the corresponding cosmologi- 
cal parameters, defining the concordance cosmological mode l 
JPage et al. 20031 ISoergel etaLlOOl ISoergel et al.'2 006). 
The universe is assumed to be globally homogeneous and 
isotropic (cosmological principle) and the CMB anisotropics 
arise from Gaussian quantum fluctuations defined in a pri- 
mordial inflationary era. From the signal processing point of 
view, the CMB is therefore a unique realization of a Gaus- 
sian and stationary (homogeneous and isotropic) random pro- 
cess on the celestial sphere and may be completely studied 
in terms of its temperature and polarization angular power 
spectra. The concordance cosmological model is defined by 
the values of the cosmological parameters giving the angu- 
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lar power spectra which best fit the experimental data. How- 
ever, among other assumptions, the basic theoretical hypothe- 
ses of Gaussianity and stationarity of the signal must imper- 
atively be thoroughly tested before the concordance model 
may be trusted. Many approaches have been proposed in 
this direction. In particular, the analysis of local features 
in the CMB temperature and polarization anisotropies might 
reveal either signatures of fundamental non-Gaussianity or 
statistical anisotropy, or foreground emissions. The scale- 
space analysis of the CMB, i.e. the combined analysis of 
scale and localization of features in the signal, is therefore 
of fundamental interest. In that perspective, the efficiency 
of the wavelet signal processing of the CMB has already 
been established. First results have also been obtained from 
the wavelet a nalysis of the WM A P all-sky temperature data, 
notably in dVielva et al. 20041 i Mukheriee & Wan g 2004t 
McE wen et al. 2005aUCruz et al. 2005ft . 

A practical approach to wavelets on the sphere was recently 
introduced in that context of scale-space analysis of the CMB 
JWiaux et al. 20051) . In the proposed formalism, the analysis 
of a signal is performed through its wavelet coefficients, re- 
sulting from the directional correlation of the signal with a 
localized filter to which dilations, rotations, and translations 
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may be applied. Analogously to the corresponding formalism 
on the plane, wavelet filters on the sphere must satisfy an ad- 
missibility condition which guarantees that any signal may be 
explicitly reconstructed from its wavelet coefficients. A cor- 
respondence principle was also established, which states that 
the inverse stereographic projection of a wavelet on the plane 
gives a wavelet on the sphere. This correspondence principle 
enables to transfer properties, such as the filter directionality 
and steerability, from the plane onto the sphere. Any non- 
axisymmetric filter is said to be directional. The good direc- 
tionality of a filter may be measured in terms of the peaked- 
ness of its auto-correlation function. The steerability of a filter 
is that property through which an arbitrary rotation of the fil- 
ter around itself may be comp uted exactl y in terms of a finite 
number of basis filters (Wia ux et al. 20051) . 

The present work develops a fast algorithm for the direc- 
tional correlation of band-limited signals and band-limited 
steerable filters on the sphere, in particular wavelet filters. 
In the context of the scale-space CMB analysis, this algo- 
rithm reduces typical computation times for the directional 
correlation of all-sky maps of several megapixels from the 
ongoing WMAP or the forthcoming Planck missions from 
years to tens of seconds, thus rendering them easily af- 
fordable. The scale-space wavelet analysis of the CMB 
temperature (T) and polarization (E and B) anisotropics, 
and in particular the identification of possible local non- 
Gaussianity or statistical anisotropy signatures, or foreground 
emissions, is thereby accessible with the same high preci- 
sion in both position and direction. The proposed algorithm 
may also be directly applied in the study of beam asymme- 
try effects on the C MB temperat ure and polar i zation dat a 
(Chall inoretal. 20001 iWandelt & Gorski 20011 N 2 2005). 
But the generic results exposed for the directional correla- 
tion on the sphere may actually find applications in many 
other fields beyond cosmology jBillow &Daniilidis 2001; 
Billow 2002). 

In §12 we define the directional correlation on the sphere as 
the scalar product of a signal and a directional filter at all pos- 
sible positions of the filter on the sphere, and for all possible 
directions at each point. It is therefore a function defined on 
the rotation group in three dimensions SO(3). By opposition, 
the standard correlation on the sphere considers a fixed direc- 
tion of the filter and is associated with a function on the sphere 
S 2 . We discuss the a priori OiL 5 ) asymptotic complexity for 
the computation of the directional correlation through simple 
quadrature, where 2L stands for the square-root of the number 
of sampling points on the sphere, also setting a band limit L 
for the signals and filters considered. The corresponding naive 
computation of the directional correlation would take several 
years on a single standard computer for fine samplings on the 
sphere of order L ~ 10 3 , and is thus inaccessible in terms of 
computation times. We then recall that the spherical harmon- 
ics and the Wigner D-functions form orthogonal bases in the 
spaces of square-integrable functions on S 2 and SO(3) respec- 
tively. We also introduce a relation giving the directional cor- 
relation as the inverse Wigner D-functions transform of the 
pointwise product between the signal and filter scalar spher- 
ical harmonics coefficients. In § [5] we recall the notion of 
steerable filters on the sphere, and give the examples of the 
first and second Gaussian derivatives. We also show that the 
directional correlation of arbitrary signals with steerable fil- 
ters is explicitly given in terms of the standard correlations 
with their basis filters. We briefly introduce spin-weighted 
spherical harmonics as orthonormal bases for the decompo- 



sition of spin functions on the sphere and derive a natural 
relation giving the standard correlation as a sum of inverse 
spin-weighted spherical harmonics transforms. In §|4] we ex- 
plicitly define our algorithm. We first discuss the global al- 
gorithmic structure, as well as pixelization choices. We then 
establish a recurrence relation which expresses spin-weighted 
spherical harmonics as linear combinations of scalar spher- 
ical harmonics. Consequently, the directional correlation of 
a scalar function with a steerable filter may be uniquely ex- 
pressed in terms of direct and inverse scalar spherical harmon- 
ics transforms. The separation of variable technique sets the 
asymptotic complexity of the scalar transform to 0(L 3 ). The 
specific use of equi-angular pixelizations allows to achieve 
the exactness of the algorithm, through a sampling theorem 
introduced by Driscoll and Healy for their fast scalar spher- 
ical harmonics transform. If the associated Legendre poly- 
nomials necessary for the transform are pre-calculated, the 
asymptotic complexity of the Driscoll and Healy transform 
on equi-angular grids drops to an optimized 0{L 2 \og^L). We 
discuss the numerical implementations of the scalar spherical 
harmonics transforms on equi-angular pixelizations (Sphar- 
monicKit package) as well as on Hierarchical Equal- Area iso- 
Latitude pixelizations (HEALPix package). A simple com- 
parative analysis of the corresponding memory requirements, 
computation times, and numerical stability shows that both 
implementations of the directional correlation algorithm are 
easily accessible on a standard computer for resolutions well 
beyond L ~ 10 3 . Finally, we briefly recall an existing alter- 
native C(L 4 ) algorithm for the directional correlation, which 
may also be optimized to an 0{L 3 ) asymptotic complexity 
through the use of steerable filters. In § |5J we consider the 
algorithm in the perspective of the scale-space analysis of the 
CMB, through the directional correlation of the temperature 
T, and of the electric E and magnetic B polarization compo- 
nents. We illustrate the performance of our algorithm through 
the explicit computation of the wavelet coefficients of a simu- 
lated CMB temperature map with the second Gaussian deriva- 
tive wavelet. A comparison of the results of the equi-angular 
and HEALPix implementations is proposed. The good pre- 
cision of the two implementations is confirmed with a clear 
advantage though for the equi-angular approach which is the- 
oretically exact. The computation times of directional corre- 
lation of maps of several megapixels are reduced from years 
to tens of seconds. We briefly conclude in §[6] 

2. DIRECTIONAL CORRELATION 
2.1. Definition 

We define here explicitly the notions of directional and stan- 
dard correlations. 

Let the function F(oj) and the filter ^(oj) be square- 
integrable functions in L 2 (S 2 ,dil) on the unit sphere S 2 . The 
point ui is given in the spherical coordinates defined in the 
right-handed Cartesian coordinate system (o,ox,oy,oz) cen- 
tered on the unit sphere as: oj = (0,ip), The angle 9 £ [0,7r] 
is the polar angle, or co-latitude. The angle <p £ [0,27r[ is 
the azimuthal angle, or longitude. The invariant measure 
on the sphere reads dVl = d cos 9d(p. Recall that an axisym- 
metric filter is by definition invariant under rotation on it- 
self. That is, when located at the North pole, an axisym- 
metric filter is defined by a function ^(9) independent of 
the azimuthal angle <p. Any non-axisymmetric filter is said 
to be directional, and is given as a general function tp) 
in L 2 (S 2 ,dil). The directional correlation is defined as the 
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scalar product between the function F and an arbitrary filter 
\P rotated on itself in a direction \ G [0,27r[, and translated 
at any point ojq = (0o,<Po) on the sphere. The Euler angles 
((po,9o,x) associated with a general rotation in three dimen- 
sions p € SO(3) represent successive rotations by \ around 
oz, do around oy, and ipo around oz. These may also be in- 
terpreted in the reverse order as successive rotations by ipo 
around oz, #o around oy', and x around oz" , where the axes 
oy' = oy'((fo) and oz" = oz"(tpo, 9q) are respectively obtained 
by the first and second rotations of the coordinate system by 
tpo and 9q (Brink & Satchler 1993). The translations of the 
filter at u> = (60, ipo) are thus associated with the two first 
Euler angles ipo around oz, 9o around oy'. The rotations of the 
filter on itself in the plane tangent to the sphere at u> = (9o, ipo) 
are rotations around oz", associated with the third Euler angle 
X- The affine transformations considered are therefore im- 
plemented through the action on the function of the oper- 
ators R(p) associated with a rotation p = (tpo, 80, x) m SO(3). 
These operators read, in the first Euler angles interpretation, 
as R(p) = RHipo)R 9 (9o)RHx) = R(^o)RHx)- The rotation 
reads [/?(/>)*](«;) = tf (fl^w). where R p = = R^ 

stands for the three-dimensional rotation matrix associated 
with p, acting on the coordinates u> = (0,ip). The result 
of the directional correlation thus explicitly gives a square- 
integrable function in L 2 (SO(3),dp) on the rotation group, for 
the invariant measure dp = dcos dodifodx- 

The directional correlation (Rty\F) of the function F by the 
filter thus explicitly reads in L 2 (SO(3),dp) as: 



(R(p)*\F) 



dn^* 



S 2 



(R- l u)F(u;). 



(1) 



The standard correlation (Rq^ \F) of F by "J, is defined by the 
scalar product between the function F and the filter "J trans- 
lated at any point u>o = (9q, tpo) on the sphere, but for a fixed 
direction, i.e. a fixed value x = 0. The result of the stan- 
dard correlation explicitly gives a square-integrable function 
in L 2 (S 2 ,dft) on the sphere: 



tf|F) = / dfitf*(tf^w)F(w). 



(2) 



The notation Rq simply identifies a three-dimensional rotation 
with x = 0. It distinguishes the standard correlation (Rq^>\F) 
from the directional correlation (R^ \F) when the arguments 
are not specified. In the particular case of an axisymmetric fil- 
ter, there is no dependence of the correlation in the filter rota- 
tion x- The directional correlation with an axisymmetric filter 
is therefore strictly equivalent to the standard correlation. 

Notice that if "J is the specific dilation by a € M* of 
a wavelet on the sphere, — > ty a , the directional corre- 
lation identifies with the wavelet coefficient of the signal, 
at the correspondin g scale: W££u iQ,x,a) = (R(pWa\F), for 
p = (<po,8o,x) (Wiauxetal. 2005). Let us finally remark 
that the correl ation is also g enerically called convolution 
on the sphere (Wandelt & Gorski 2001). We use a specific 
terminolog y to avoid any confusion with other definitions 
JPriscoll & Healv 1994tlHealv et al. 20031 

2.2. A priori computation cost 

The numerical computation cost associated with a naive 
discretization, or quadrature, of the relation Q for the direc- 
tional correlation of functions on the sphere may easily be 
estimated in the following way. 



First, let N p = (2L) 2 represent the number of sampling points 
ui in a given pixelization of the sphere. The quantity 2L repre- 
sents the mean number of sampling points in the position vari- 
ables 9 and tp. A simple extrapolation of the Nyquist-Shannon 
theorem on the line intuitively associates L with the band 
limit, or maximum frequency, accessible in the "Fourier" in- 
dices conjugate to 9 and cp for the signals and filters consid- 
ered. To an equivalent sampling on 2L points in the direction 
X also corresponds a band limit L in the conjugate Fourier 
index. 

Second, considering a simple discretization of the relation 
of directional correlation, each integration on to = (9, tp) on the 
sphere, a scalar product, has an asymptotic complexity 0(L 2 ). 
The directional correlation is associated with the calculation 
of a scalar product for each discrete p = (ipo,9o,x) on SO(3). 
The corresponding naive asymptotic complexity is therefore 
of order 0(L 5 ). 

Third, we consider fine samplings of several megapixels on 
the sphere associated with a band limit around L~ 10 3 . In par- 
ticular, present CMB experiments such as the WMAP mission 
provide all-sky maps of around three megapixels. For such a 
fine sampling, the typical computation times for (2L) 2 mul- 
tiplications and (2L) 2 additions are of order of 0.03 seconds 
on a standard 2.2 GHz Intel Pentium Xeon CPU, for double- 
precision numbers. We take this value as a fair estimation 
of the computation time required for a scalar product in Q. 
Consequently, a unique 0(L , ) directional correlation would 
take several years at that band limit L ~ 10 3 on a single stan- 
dard computer. Moreover, depending on the application, the 
directional correlation of multiple signals might be required. 
Typically, thousands of simulated signals are to be considered 
for a statistical (CMB) analysis. 

In conclusion, the computation time for a simple 0(L ) dis- 
cretization of the relation for the directional correlation of 
functions on the sphere is absolutely unaffordable for fine 
samplings with a band limit around L ~ 10 3 in 9, tp, and %. 
This conclusion remains when the use of multiple computers 
is envisaged. It is even strongly reinforced in the perspec- 
tive of a scale-space analysis on the sphere from finer pix- 
elizations. In particular, the future all-sky CMB maps of the 
Planck mission will reach samplings of fifty megapixels, i.e. 
I. -4 • 10 ! . 

2.3. Directional correlation as a Wigner transform 

We here discuss the decomposition of the directional corre- 
lation as an inverse Wigner D-functions transform on 5(9(3). 

First, we recall the decomposition of functions on S 2 and 
5(9(3) in standard scalar spherical harmonics and Wigner D- 
functions respectively. The standard scalar spherical harmon- 
ics Yi m (u>), with I £ N, m G Z, and \m\ < I, form an orthonor- 
mal basis for the decomposition of func tions in L 2 (S 2 , dil) on 
the sphere, with ui = (0,ip) ( Varshalovich et al. 1989). They 
are explicitly given in a factorized form in terms of the asso- 
ciated Legendre polynomials P'"{cos 9) and the complex ex- 
as 



ponentials e' mip 

Ylm(9,p) = 



21+1 (l-m)\ 



1/2 



Pf (cos 9) e" 



(3) 



47r (l + m)\ 

This corresponds to the choice of Condon-Shortley phase 
(-1)"' for the spherical harmonics, ensuring 



for 

(-i) m Y; m (to) 



the relation 
K;(_ m) (w). This phase is here included 



in the definitio n of the associated Le gendre polynomials 
( Abramowi tz & Stegun 19651 IVarshalovich et al. 19891) . An- 
other convention JBrink & Satchler 19931) explicitly transfers 
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it to the spherical harmonics. The orthonormality and com- 
pleteness relations respectively read: 



.s- 



dQ Y* m (uj) Y Vm i (lo) = 6u> 5„ 



and 



^^Yf m {L,')Y lm (u,) = 5{u,'- 

IE® \m\<l 



(4) 



(5) 



with 5(lo' -uo) = S(cos 9' -cos 6)8(tp' - tp). Any function G{uf) 
on the sphere is thus uniquely given as a linear combination 
of scalar spherical harmonics (inverse transform): 



G(UJ) = ^ Gl m Yl m (UJ) , 

/eN \m\<l 



(6) 



for the scalar spherical harmonics coefficients (direct trans- 
form) 

G lm = f dVLY* m (uj)G(uj), (7) 

with \m\ < I. 

The Wigner ^-functions D l mn (p), for p = (tp, 6, \) € SO(3), 
and with / e N, m,n G Z, and |w|,|n| < I, are the ma- 
trix elements of the irreducible unitary representations of 
weight I of the rotation group SO(3), in L 2 (SO(3),dp). By 
the Peter- Weyl theorem on compact groups, the matrix el- 
ements D l *„ also form an orthogonal basis in L 2 {SO(3),dp) 
dVarshalovich et al. 19 89). They are explicitly given in a fac- 
torized form in terms of the real Wigner li-functions £/,'„,(#) 
and the complex exponentials, e~"" v and e~'" x , as 

D l mn (vAx) = e- im *d l mn (6)e- in x. 
The Wigner ^-functions read 



(8) 



Al ,m ^(-l)'[(l + my.(l-m)\(l + n)l(l-ny.]^ 2 

t=Ci 



{l + m-t)\{l-n-t)\t\(t + n-m)\ 



n ,-\2l+m-n-2t , . „ , v 2(+«-m 

iO/2) (smv/2) , 



(9) 



with the summation bounds C\ = max(0,m-n) and C2 = 
min{l + m,l-n) defined to consider only factorials of posi- 
tive integers. They satisfy various sym metry properties on 
their indices ( Varshalovich et al. 1989). The orthogonality 
and completeness relations of the Wigner D-functions respec- 
tively read: 



dpD' mn ( P )D' m r„,(p): 



and 



50(3) 



Ei^ E D' mn (p')D' l : n (p)=s(p'-p), (id 



21+1 



(10) 



let 



\m\,\n\<l 



with 5(p' — p) = S(ip' - tp)S(cos 9' - cos 9)5(x' - X)- Any func- 
tion G(p) in L 2 (SO{3),dp) is thus uniquely given as a linear 
combination of Wigner D-functions (inverse transform): 



^ 8tt 2 



(12) 



\n\<l 



for the Wigner D-functions coefficients (direct transform) 



dpD' mn (p)G(p), 



(13) 



with \m\, \n\ < I. 

Second, we expose the decomposition of the direc- 
tional correlation in Wigner D-functions. The Wigner D- 

functions coefficients (Rty\F) nm of the directional correlation 
(R(p)^>\F) are given as the pointwise product of the scalar 

spherical harmonics coefficients F/ m and *&* n . The following 
directional correlation relation indeed holds: 



(R(p)*\F)=^ 



/en 



2/+1 
8^2" 



E 

n I . I n I < / 



{R^\F) mn D l ; m {p), (14) 



with the Wigner D-functions coefficients on SO(3) (R*f?\F 
given as 

- ' 8?r 2 - - 



mn m „ = 



21+1 



(15) 



The derivation of this result goes as follows. The orthonor- 
mality of scalar spherical harmonics implies the following 

Plancherel relation (R*|F) = E/eN E| m |</^L^»- The ac- 
tion of the operator R(p) on a function G(uj) in L?{S 2 ,dVl) 
on the sphere reads in terms of its scalar spherical harmon- 
ics coefficients and the Wigner D-functions as: [R(p)G]i m = 
E| n |</ D l m „(p)Gi„. Inserting this last relation for ^ in the for- 
mer Plancherel relation finally gives the result. In the par- 
ticular case of an axisymmetric filter ^(9), the above rela- 
tion reduces to the following standard correlation relation: 

(Ro^\F) lm = 2ir i &*Fi m , where <3>/ stands for the Legendre co- 
efficient of the filter. 

3. FILTER STEERABILITY 
3.1. Definition 

We here recall the notion of filter s teerabili ty on 
the sphere iWiaux et al. 2005t iFreeman & Adelson 199 It 
iSimoncelli et al. 19921) and give explicit examples. A direc- 
tional filter iff in L 2 (S 2 ,dil) on the sphere is steerable if any 
rotation by \ G [0,27r[ of the filter around itself R : (x)^ may 
be expressed as a linear combination of a finite number of ba- 
sis filters ^ff m : 



M 



(16) 



m- 1 



50(3) 



The weights k m (x), with 1 < m < M, and M £ N, are called 
interpolation functions. In particular cases, the basis filters 
may be specific rotations by angles \m of the original fil- 
ter: ^f m = R z (xm)^- Steerable filters have a non-zero angular 
width in the azimuthal angle tp which makes them sensitive 
to a whole range of directions and enables them to satisfy the 
relation dl6> . In the spherical harmonics space, this non-zero 
angular width corresponds to an azimuthal angular band limit 
N in the index n associated with the azimuthal variable tp: 

$/„ =0 for \n\>N. 

The inverse stereographic projection on the sphere of the 
N'J 1 derivative in direction x of radial functions on the plane 
(tangent at the North pole) are steerable filters. They may 
be rotated in terms of M = Nd + 1 basis filters, and are band- 
limited in tp at N = Nd + 1 ■ We give here the explicit examples 
of the normalized first and second Gaussian derivatives. A 
first derivative has a band limit N = 2, and only contains the 
indices n = {±1}. It may be rotated in terms of two specific 
rotations at x = and x = 7r /2, corresponding to the inverse 
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projection of the first derivatives in directions i and y, "J °* and 
ty d f respectively: 

[rHx)® 9 *] (w) = **(w)cosx+* 9t Msinx- (17) 

The normalized first derivatives of a Gaussian in directions i 
and y read: 



- I 1 - (an 2 r ) e - 2tan2(e / 2) tan - cos y> 



$ 9 *™ ! » 0,<p)=Jl(l+ tan 2 e - 2t - 2 ( / 2 ) tan | 



sin ip. 
(18) 



A second derivative has a band limit N = 3, and contains the 
frequencies n = {0,±2}. It may be rotated in terms of three 
basis filters. It reads indeed in terms of the inverse projection 

ft 2 r) 1 

of the second derivatives in directions i and y, w * and f j 
respectively, and the cross derivative ty®* 8 ) as: 



R Z (X)^ C 



(uj)=^ 9 ' (w)cos 2 x+ (w)sin z x 
+^r% a f( w ) sin2x- 



(19) 



The correctly normalized second derivatives of a Gaussian in 
directions x and y read: 



d~ (gauss) 




— I 1+tan 2 ^ 1 e 



-2tan 2 (fl/2) 



# a j(s fl " si >(0,^) = 




1 -4 tan 2 — cos 2 ip 



, 1+tan 2 - e - 2tan(e/2) 
3tt V 2. 



l-4tan 2 -sin 2 (y5 



1+tan 2 - e " 2tan2(e/2) 
/3tt V 2,' 



tan - sin2(/) 
2 r 



3.2. From directional to standard correlation 



(20) 



We here express the directional correlation of arbitrary sig- 
nals with steerable filters in terms of standard correlations. 

The relation of steerability d!6i is obviously preserved by 
linear filtering. The directional correlation of a signal F by a 
steerable filter W therefore satisfies the same steerability rela- 
tion as the filter itself: 



{R (p) *\F) = k >» (X) (fi (^o) * ffl \F) , 

m= 1 



(21) 



for p = (ipo,8o,X)> an d w o = (<Po,@o)- The steerability therefore 
enables the computation of the directional correlation of a sig- 
nal F with a steerable filter as a linear combination with 
known weights k m (x), of M standard correlations with the ba- 
sis filters ^,„. For steerable filters with M <<L, this quantity 
M disappears from asymptotic complexities calculations. In 
conclusion, the asymptotic complexity of the directional cor- 
relation with a steerable filter is reduced to the asymptotic 
complexity of a standard correlation, to which must be added 
the 0(L?) operations required for the explicit computation of 



the directional correlation through the relation \2\\ . Already 
notice that the a priori computation cost associated with the 
discretization of the directional correlation integral Q with 
a steerable filter is lowered from 0{L 5 ) to 0(L 4 ) through the 
calculation of standard correlation integrals (|2j with the cor- 
responding basis filters. 

3.3. Standard correlation as a spin-weighted spherical 
harmonics transform 

We here study the decomposition of the standard correla- 
tion as a sum of inverse spin-weighted spherical harmonics 
transforms on the sphere. 

Let us first recall the notion of spin n function in L 2 (S 2 ,dVl) 
on the sphere and the related decomposition in spin-weighted 
spherical harmonics of spin n. As defined in § |2] in the co- 
ordinate frame (o, ox, oy, oz), the Euler angles (tp,9,x) associ- 
ated with a general rotation in three dimensions may be in- 
terpreted as successive rotations by tp £ [0,27r[ around oz, 
9 £ [0,7r] around oy'(ip), and \ <= [0>2tt[ around oz"(<p,9). 
The local rotations of the basis vectors in the plane tangent 
to the sphere at oj = (9, cp) are rotations around oz", asso- 
ciated with the third Euler angle \. Spin n functions on 
the sphere „G(uj), with n £ Z, are defined relatively to their 
behaviour und er the cor responding right-handed rotations 
by Yn as (|Newman & Penrose 1966; Goldb erg et al. 19671 
Carm eli 19691) : 

(22) 



„G'(w) = ^%G(w). 



We emphasize that the rotations considered are local trans- 
formations on the sphere around the axis oz" = oz"(<p,0), af- 
fecting the coordinate x in the tangent plane independently 
at each point oj = (6,<p), according to x' = X~Xo- They 
are to be clearly distinguished from the global rotations R z x 
associated with the alternative Euler angles interpretation, 
which affect the coordinates of the points oj = (9, ip) on the 
sphere. Our sign convention in the exponential is coher- 
ent with the definition (I23> here below for the spin-weighted 
spherical harm onics of spin n. It is op posite to the orig- 
inal definition ( Newman & Penrose 1966), while equivalent 
to r ecent notations used in the context of the CM B analy- 
sis dZaldarriaga & Seljak 1997; Challi noret al. 2 000). Re- 
calling the factorized form (|8j, spin functions are equiva- 
lent^ defined as the evaluation at x = of any function 
in L 2 (SO(3),dp) resulting from an expansion for fixed in- 
dex n in the Wigner D-functions D' mn (ip,8,x)- The functions 
D l mn (ip,9,0) or D l * { _ n) (ip,8,0) thus naturally define for each n 
an orthogonal basis for the expansion of spin n functions in 
L 2 (S 2 ,dfl) on the sphere. Their normalization in L 2 (S 2 ,dfl) 
defines the spin-weighted spherical harmonics of spin n: 



n Y lm (9,<p) = (-iy 



21+1 



4n 



-D , : { _ n) (ip,9,Q), 



(23) 



with I £ N, / > \n\, and m £ Z, \m\ < I. They are thus ex- 
plicitly given in a factorized form in terms of the real Wigner 
(i-functions d' mn (9) and the complex exponentials e" nv as 



n Y lm {9,p) = (-!)" 



21+1 
4n 



d'm(-n) (0) 



3 im<p 



(24) 



In particular, the symmetry properties of the Wigner d- 
functions imply the generalized relation (— l) n+m „Y^(u/) = 
_„7;(_ m )(w). The spin spherical harmonics explicitly iden- 
tify with the standard scalar spherical harmonics: oYi m (ui) = 
Y lm (oj), through the relation d\jjS) = [ l #J?(cos 9). The 
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orthonormality and completeness relations respectively read 
from relations (II Oi and il It . as 



S 2 



dQ n Y* m {uf) n Yl'm' (W) = Sll'S„ 



(25) 



and 



" Y '>" ( w » Y "» (w) = 5 ( w ' - w ) ' (26) 

;SN |m|</ 

with <5(cj' — uS) = i5(cos 0' - cos 9)8(ip' - p). Any spin n function 
„G(w) on the sphere is thus uniquely given as a linear combi- 
nation of spin n spherical harmonics (inverse transform): 



„G(w) = E "^' m nYlm ( w ) > 
/GN |m|</ 



(27) 



for the spin-weighted spherical harmonics coefficients (direct 
transform) 

/ <ffi„Y,* (w)G(w), (28) 

with I > \n\, and \m\ < I. 

The directional correlation at \ = identifies with the stan- 
dard correlation. From the above definitions, the relation dl4> 
therefore leads to the following standard correlation relation 
on the sphere, expressed as a sum of inverse spin-weighted 
spherical harmonics transforms: 



E (Ro^\F) lmnn Y lm (^) 



t> n \m\<l 



, (29) 



with coefficients (Ra^lF) lmn explicitly defined as the point- 
wise products 



» = (-D" 



4-TT 

21+1 



(30) 



For each n, these coefficients may be understood as spin- 
weighted spherical harmonics coefficients, but to be clearly 
distinguished from the spin-weighted spherical harmonics co- 
efficients „(R(^\F) lm of the standard correlation, which do 



not find a pointwise product expression. Considering a real 
filter also implies the further simplification: (— 



/(-«) ' 



4. FAST ALGORITHM 



4. 1. Algorithm and pixelization 

In the following, we define the global structure of our algo- 
rithm and discuss pixelization choices. 

In terms of the relations (I29> and d30i . the computation of 
the standard correlation of a signal F by a steerable filter ^ 
may be performed with the following global structure: a di- 
rect scalar spherical harmonics transform of the signal and 
the filter, and F/ m , a correlation in spherical harmonics 
space through a pointwise product, and a sum of inverse spin- 
weighted spherical harmonics transforms to obtain (R^\F). 

We consider band-limited signals F at some band limit L 
on the sphere S 2 , that is by definition Fi m = for I > L. From 
( I30> . the standard correlation of a band-limited signal F with 
a band limit L by an band-limited filter \fr on the sphere is thus 



also band-limited, at the same band limit: (Rn^lF) 



:0for 



M)*!* //to,- 

I > L. In order to achieve the limit frequency L, an extrapo- 
lation of the Nyquist-Shannon theorem on the line typically 
requires pixelizations with at least 2L sampling points in the 



polar angle 6 S [0,7r] and the azimuthal angle tp G [0,27r[ on 
the sphere. In the context of the signal processing of the 
CMB maps, many pixelization schemes have been consid- 
ered. We only quote here the HEALPix scheme ( Hierarchi- 
cal Equal Area iso-Latitude Pixelization) (G orski et al. 200.51 
notably used for the WMAP and Planck experiments, and 
the GLESP pixelization (Gauss-Legendre Sky P ixelization) 
JPoroshkevich e t al. 20053 iDoroshkevich et al. 2005bJ. In 
particular, on a 2L x 2L equi-angular grid, a sampling re- 
sult on the sphere states that the direct scalar and spin- 
weighted spherical harmonics coefficients of a band-limited 
function on the sphere may be computed exactly up to a 
band limit L as a finite weighted su m, i.e. a quadrature, of 
the sampled values of that function (Drisco U~& Healv 19941 
iKostelec et al. 20 001 The weights are defined from the struc- 
ture of the Legendre polynomials P/(cos60 on [0,7r]. They 
are functions of 9, but are independent of the azimuthal an- 
gle ip and of the spin n considered. The theoretical exactness 
of computation represents an advantage of equi-angular grids 
relative to other samplings. Inverse scalar and spin-weighted 
spherical harmonics transforms of band-limited functions can 
obviously be evaluated exactly on any grid as they are explic- 
itly defined as finite sums from relations (|6j and J27i with 
1<L. 

4.2. From spin-weighted to scalar spherical harmonics 

We propose here an original expression of spin-weighted 
spherical harmonics as simple linear combinations of scalar 
spherical harmonics. 

The spin n spherical harmonics may be related to spin 
n ± 1 spherical harmonics thro ugh the action of spin 
raising and lowering operators (Newman & Penrose 1966; 
Goldberg e t al. 196711 . The action of the spin raising <3 and 
lowering 3 operators on a spin n function „G, giving a spin 
n + 1 and n - 1 function respectively, is defined as 



[S„G](^): 



and 



[5 n G](0,<p) = 







d 



89 sin 9 dip 



-sin" ft | — + - — - — ] sm 9 „G (6,<p), 

(31) 



d_ 

89' 



8 



sin 9 8<p 



sin" 6 n G 



(0,P), 
(32) 

with, under rotation by xo m the tangent plane at u> = 
(9,p): [i3„G]V> = e- i( " +1 >x°[9„G]M and [8„G]'M = 
e ,( " 1)xo [3 „G](lj). In these terms, the spin-weighted spher- 
ical harmonics of spin n are related to spin-weighted spher- 
ical harmonics of spin n + 1 and n - 1 through the following 
relations: 



[3 „y /m ]M = [(/-«)(/+«+!)] 



1/2 



Y\Ylm (w), 



(33) 



and 



[3 n Y lm ] (lj) = -[(/ + «)(/-»+ 1)] 1/2 n -iY lm (u) , 

also implying 

[33 n Y lm ] (uj) = -(l-n)(l + n+l) n Y lm (to) . 

The corresponding direct relation between the spin-weighted 
spherical harmonics of spin n and scalar spherical harmonics 
reads: 



(34) 



(35) 



(/-»)! 
(l+ri)\ 



1/2 



[B"Y lm ] QJ), 



(36) 
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for < n < I, and 

tiYlm (w) = 



(/ + «)! 
(/-«)! 



1/2 



(-D" rr-yJcw), 



(37) 



for-/<n<0. 

Beyond these standard relations, we establish a recurrence 
relation free of derivatives for spin-weighted spherical har- 
monics as follows. Derivative relations s atisfied by the 
Wigner D-functions ( Varshalovich et al. 1989) allow to ex- 
change the derivative in 9 in J33i and J34i for simple linear 
combinations. This gives the spin n spherical harmonics „Y/ m 
as a linear combination of the spin n =p 1 spherical harmonics 
n^fiYim and n zpiT(/-i) m : 



n Yi m (9,f) = 



1/2 



m^flcosO 



l±n J Isint 

(l±n-l)(l 2 -n 
T±~n 



21+1 
21-1 



1/2 



Zsin6> 



n^lY(!-\) m 



(38) 

under the convention that „Yi„ t is defined to be zero for / < 
max(|m|, \n\). Grouping terms through their l/sinf? or cot 9 
pre-factors on the sphere, one gets in compact form the fol- 
lowing 2-terms relation: 



n Y, m (9,ip) = 



T a (l, 1 ) ± P(l. 



mn) 



sine 



■Ta ( t)C0tf 



„ T iY lm (0^), 



(39) 

with rv ± - /- 'Tn+U l/2 anH fl± _ l f M+l g±«-l)g 2 -m 2 ) -.l/2 

with a (/n) - (-^-) , and p ((mj!) - ^ J ■ 

The operators S ±l : I — > / ± 1 define a one-unit shift in 
the index Z: „ T iF/ m ](w) = „ T iT(/-i) m (w). By a sim- 
ple |«|-steps recurrence, spin-weighted spherical harmonics 
may therefore by expressed as a 2l"l-terms linear combi- 
nation of scalar spherical harmonics. The same relation 
between spin-weighted and scalar spherical harmonics may 
also be obtained from relations (I36i and (I37> and direct 
derivative relations for the as sociated Legendre polynomials 
jAbramowitz & Steeun 19651) . 

4.3. Detailed algorithmic structure from scalar spherical 
harmonics transforms 

In the following, we discuss the detailed structure of the 
algorithm, which may essentially be decomposed in a combi- 
nation of scalar spherical harmonics transforms, and study the 
corresponding asymptotic complexity. 

From the recurrence (I39i . the inverse spin-weighted trans- 
forms of (Rq^\F) , at each n required in the relation \29\ 
for the standard correlation may be decomposed as a linear 
combination of inverse scalar spherical harmonics transforms. 
Hence, the complete standard correlation of a band-limited 
signal with a band-limited steerable filter essentially relies on 
direct and inverse scalar spherical harmonics transforms. 

The a priori complexity associated with the naive com- 
putation of the direct scalar transform integral (Q in (9,ip) 
on the sphere through simple quadrature, for all (l,m) with 
\m\ < I < L, or with the sum on (/,m) for all discrete points 
(9,ip) in the inverse scalar transform (|6j is naturally of or- 
der 0(L 4 ). The separation of variables (0 in the scalar 
spherical harmonics into the associated Legendre polynomials 
P™(cos0) and the complex exponentials e' mip allows to com- 
pute direct and inverse scalar spherical harmonics transforms tation and algorithm implementation (SpharmonicKit package). 



as successive transforms in each of the variables 9 and ip in 
0(L?) operations. The application of the concept of separa- 
tion of variables on the sphere only requires that the sampling 
in the polar angle 9 be independent of the azimuthal angle 
ip. This criterion is met for a large variety of pixelizations. 
We will consider the corresponding sta ble numerical imple- 
mentations on HEALPix pixelizations (Gorskietal. 2005 ) 2 . 
For a 2L x 2L equi-angular grid in {9, if) on the sphere, a fur- 
ther op timized algorithm exist s which is due to Driscoll and 
Healy (Driscoll & Healv 1994). It explicitly takes advantage 
of a recurrence relation in I on the associated Legendre poly- 
nomials P™(cos 9) to compute the associated Legendre trans- 
forms in ©(Llogj L) operations for each m. The Fourier trans- 
forms in e' mip are computed in C(Llog 2 L) operations for each 
9 through standard Cooley-Tukey fast Fourier transforms. In 
these terms, the direct and inverse scalar spherical harmon- 
ics transforms of a band-limited function of band limit L on 
the sphere are computed in 0(L 2 log 2 L) operations. More- 
over, contrarily to the HEALPix implementation, the algo- 
rithm is theoretically exact. The exactness of the direct trans- 
form relies on the specific choice of weighting functions only 
depending on 9, as required by the sampling result on equi- 
angular grids on the sphere. Again, the exactness of the eval- 
uation of the inverse scalar spherical harmonics transform of 
a band-limited function relies on the fact that is expressed as 
the simple finite sum (|6) with / < L. We will consider the cor- 
responding stable numerical implementations of the Sphar- 
monicKit package ( Healv et al. 20031: iHealv et al. 2004F . As 
detailed in § 14.41 if high resolutions are to be considered 
(L > 1024), the associated Legendre polynomials are not 
pre-calculate in order to avoid a too large RAM memory 
consumption. They must be computed on the fly and the 
corresponding 0(L y ) time consumption is included in the 
global asymptotic complexity for the spherical harmonics 
transform. Both the C(L 3 ) HEALPix implementation and 
the 0(L 2 logj L) SpharmonicKit implementation therefore ef- 
fectively have the same global 0(L y ) asymptotic complexity 
when the calculation of the associated Legendre polynomials 
is considered. 

Each inverse spin-weighted transform required in (I29i a 
priori reads as 2'"l-terms linear combination from the recur- 
rence ( 1391 . But these terms may be grouped in \n \ + 1 in- 
verse scalar spherical harmonics transforms with pre-factors 
cot 1 ' 9/ sin 9 (9 for p,q £ N and p + q= \n\. In addition, for each 
p and q with p + q= \n\, the coefficients for the spins n and — n 
may be grouped before applying the corresponding inverse 
scalar spherical harmonics transform. This allows to consider 
only spin-weighted transforms for positive spins and corre- 
spondingly further reduces the number of scalar transforms 
required for the standard correlation. 

Notice that another recurrence relation might be substituted 
to ( I39l > for the treatment of the inverse spin-weighted trans- 
forms through inverse scalar transforms. This other relation 
explicitly relates „Y lm with „^Y hn , „^Yq-u m , a nd „^\YQ+i )m 
( Varshalovi ch et al. 1989t iKostelec et al. 20001) . The term 
n+-iY(i+i)m(u) = [S +l n^iYim](u) actually raises the band limit 
of the associated scalar functions to be analyzed to L+ \n\ 
after the \n\ -steps recurrence leading from spin- weighted to 
scalar spherical harmonics, with \n\ < L. From a practical 



2 See healpix.jpl.nasa.gov for documentation and algorithm implementa- 
tion (HEALPix2. 1 software). 

3 See www.cs.dartmouth.edu/~geelong/sphere/ for unpublished documen- 



8 



point of view, in a HEALPix scheme for a fixed resolution, 
this will reduce the numerical precision of computation, as 
the implementation is known to make larger errors at higher 
I's, On 2L x 2L equi-angular grids, the SpharmonicKit pack- 
age is technically limited to consider coefficients up to L, and 
numerical errors will occur due to the absence of considera- 
tion of the coefficients between L and L+ \n\. No such issue 
occurs from relation d39l which preserves the band limit L to 
be considered for the associated scalar functions. 

In that context, let us analyze the overall asymptotic com- 
plexity of the standard correlation algorithm proposed. We 
consider steerable filters with low azimuthal angular band 
limit N « L, which constrain the maximum spin value to 
\n\ < N. Typically N = 2 for a first Gaussian derivative, and 
N = 3 for a second Gaussian derivative. The index n there- 
fore decouples from any asymptotic complexity contribution. 
The required direct and inverse scalar transforms set the over- 
all asymptotic complexity of the standard correlation to 0(L?) 
(and even to 0{L 2 \og\L) on equi-angular grids if the associ- 
ated Legendre polynomials were to be pre-calculated). The 
asymptotic complexities of additional contributions to the al- 
gorithm are of order OiL 2 ) and therefore negligible. First, 
the pointwise product d30i required after the direct transforms 
of the signal and the filter is indeed of order Oil}). Sec- 
ond, the grouping of the 2'"' terms in the recurrence for each 
inverse spin-weighted spherical harmonics transform at spin 
n, as \n\ + 1 terms on which inverse scalar spherical harmon- 
ics are computed, contribute to the asymptotic complexity as 
2'"' x 0{L 2 ). This contribution is also negligible for a first or 
second Gaussian derivative filter. 

In summary, our algorithm for the standard correlation of 
a band-limited signal with a band-limited steerable filter is 
structured as follows, with related asymptotic complexities: 

• Direct scalar spherical harmonics transforms, and 
F[ m : C(L 3 ) on all pixelizations (0(L 2 log 2 , L) on an equi- 
angular grid if the associated Legendre polynomials are 
pre-calculated). 

• Correlation (Rq^\F) lmn in spherical harmonics space: 
0(L 2 ) through (HJ. 

• Sum of inverse scalar spherical harmonics transforms 
to obtain (R ^\F) through l|29) and l|39): 0(L 3 ) 
on all pixelizations (0(L 2 log 2 L) on an equi-angular 
grid if the associated Legendre polynomials are pre- 
calculated). 

Let us finally notice that the reduction of complexity to 
0(L 2 log 2 L) thanks to the use of an equi-angular grid beyond 
the calculation of the associated Legendre polynomials pro- 
vides anyway an asymptotic reduction of the overall compu- 
tation times for the corresponding correlation. It is therefore 
already of interest only in that respect. But the explicit cal- 
culation of the directional correlation from the relation \2\\ 
in terms of the computed standard correlations requires an 
additional OiL 1 ) operation. However, we emphasize that all 
the information for the directional correlation of a signal with 
a steerable filter is already contained in the standard corre- 
lations with its basis filters. In many applications the ex- 
plicit computation of the directional correlation through the 
relation ( 12 II is not required. The interest may reside in the 
computation of the values of the directional correlation in a 



small number of specific directions Xi at each point luq on the 
sphere. One may also want to determine the direction x at 
each point cjq whic h maximizes th e value of the directional 
correlation. In (Wiaux et al. 2005), this last application is il- 
lustrated in the perspective of the detection of the precise di- 
rection of local features in the CMB temperature and polariza- 
tion anisotropics, through a steerable wavelet analysis. These 
features may represent potential signatures of fundamental 
non-Gaussianity or statistical anisotropy, or foreground emis- 
sions. The overall computation times for such a directional 
analysis are only driven by the asymptotic complexity of the 
standard correlations. 

4.4. Equi-angular and HEALPix implementations 

We concisely compare the numerical SpharmonicKit and 
HEALPix packages for the scalar spherical harmonics trans- 
forms, in terms of memory requirements, computations times, 
and numerical stability. Those properties of the scalar spher- 
ical harmonics transforms implementations characterize the 
overall standard correlation implementations themselves. We 
also comment on the issue of the change of the pixelization 
on which the data are sampled, considering equi-angular and 
HEALPix grids. 

Calculations are performed on a 2.20 GHz Intel Pen- 
tium Xeon CPU with 2 Gb of RAM memory. Ran- 
dom real test-signals are considered, with band limits L e 
{128,256,512,1024}. Without loss of generality, these real 
test-signals are defined through their scalar spherical har- 
monics coefficients G/ m with \m\ < I < L , with independent 
real and imaginary parts uniformly distributed in the interval 
[-1,+1]. The reality condition G/ ( _ m) = (-l) m G* m obviously 
reduces computation times by a factor two relative to generic 
complex signals. The inverse and direct scalar spherical har- 
monics transforms are successively computed, giving numer- 
ical coefficients Hi m . The SpharmonicKit implementation is 
coded in C. The corresponding inverse transform is computed 
on 2L x 2L equi-angular pixelizations. For L = 1024, this 
corresponds to maps with N P i X = 4194304 pixels. The direct 
transform is then recomputed exactly, i.e. to the computer's 
numerical precision, up to the band limit L. The HEALPix 
implementation considered is coded in Fortran90. The cor- 
responding inverse transform is computed with a resolution 
identified by the parameter N S id e = L/2, identifying maps with 
Npix = \2N 2 ide equal-area pixels. For L = 1024, the value 
Nside = 512 defines maps with N P i X = 3145728 pixels. The di- 
rect transform is then recomputed with a very good accuracy 

Up tO L = 2N s ide- 

RAM memory requirements for the scalar transforms are 
as follows. If the required associated Legendre polynomi- 
als Pj"(cos9) are pre-calculated once for all values of /, 9, 
and m, and stored in RAM memory, the number of real val- 
ues of associated Legendre polynomials P"'(cos9) stored in 
RAM memory for all I, 9, and m is of order O0). The 
overall memory requirements allowing both direct and inverse 
transforms are still easily accessible but already higher than 
1 Gb for L = 1024, for double-precision numbers. The pre- 
calculation computation time itself is of order 0(L 3 ) through 
the use of a recurrence relation in I on the associated Leg- 
endre polynomials. As a pre-calculation, the correspond- 
ing time consumption is not to be taken into account in the 
reported computation times, which consequently remain of 
order 0(L 2 \og\L) and 0(L 3 ) in the SpharmonicKit and the 
HEALPix implementations respectively. However, the corre- 
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sponding memory requirements rapidly become unaffordable 
for resolutions higher than L = 1024, such as the resolution 
expected for the Planck experiment. We therefore chose the 
SpharmonicKit and HEALPix implementations for which the 
required associated Legendre polynomials Pj"(co&9) are cal- 
culated on the fly, and stored in RAM memory for all val- 
ues of I and 9, but for only one value of m at the time. The 
reported computation times therefore inevitably include the 
0(L?) calculation of the Legendre polynomials, to be added to 
the remaining 0(L 2 log 2 L) and 0(L?) terms in the Spharmon- 
icKit and the HEALPix implementations respectively. Both 
implementations effectively acquire the same global 0(L?) 
asymptotic complexity. For fixed m, the number of real val- 
ues P ; '"(cosfT) stored in RAM memory for all values of I and 
9 is 0(L 2 ). For the SpharmonicKit and HEALPix implemen- 
tations, the overall memory requirements allowing both direct 
and inverse transforms are respectively of order 1 .3 x 10 2 Mb 
and 4.6 x 10 1 Mb for L = 1024, for double-precision num- 
bers. They are therefore significantly lower than the corre- 
sponding values when the associated Legendre polynomials 
are pre-calculated, and Planck-type resolutions are easily af- 
fordable. 

The related computation times associated with these global 
0(L?) SpharmonicKit and HEALPix implementations for the 
inverse and direct transforms are given in Table They 
are averages over 5 random band-limited real test-signals. 
They are of the order of tens of seconds at L = 1024 both 
for the SpharmonicKit and the HEALPix implementations, in 
rough agreement with our previous intuitive estimation that 
an 0(L ) operation requires 0.03 seconds on our standard 
computer. For the SpharmonicKit implementation, compu- 
tation times evolve from 1.0 x 10 _1 seconds for L = 128 to 
4.6 x 10 1 seconds for L = 1024. For the HEALPix imple- 
mentation, notice that an iterative scheme can be used for 
the direct transform, in order to enhance the accuracy of the 
calculation. After several iterations, the process converges 
towards a limit precision, which however never reaches the 
precision of the SpharmonicKit approach. For zero iteration 
(i = 0), computation times evolve from 4.8 x 10~ 2 seconds for 
L = 128 to 1.5 x 10 1 seconds for L = 1024. The HEALPix 
implementation with a unique iteration is therefore slightly 
more rapid than the SpharmonicKit one, at least up to the 
band limit L = 1024. Computation times grow significantly 
with the number i of iterations used. For one iteration (i =1), 
the computation times for the direct transform are already big- 
ger than for the SpharmonicKit implementation. They evolve 
from 2.7 x 10"' seconds for L = 128 to 7.8 x 10 1 seconds for 
L= 1024. 

The relative maximum error on the coefficients associ- 
ated with the scalar spherical harmonics transforms is de- 
fined as max/ ,,, |(G/,„ - i/; m )/G/ m |, where | • | here denotes 
the complex norm. The relative root mean square numer- 
ical error is defined as the ratio of the L 2 -norm of G — H 
and G, that is through the Plane herel relation: (Xw m \Gi m — 

Him\ 2 /Xwjn \Gim\ 2 )^ 2 - The numerical errors given in Table|2] 
are averages for transforms over 5 random band-limited real 
test-signals. For the SpharmonicKit implementation, errors 
are extremely small as expected from the theoretical exactness 
of the algorithm. Relative maximum errors on the coefficients 
and relative root mean square errors are indeed bounded by 
1 .8 x 10" 7 and 7.5 x 10" 10 respectively, up to L = 1024. For the 
HEALPix implementation with zero iteration (i = 0), the rela- 
tive maximum errors on the coefficients can reach 2.9 x 10 1 , 



TABLE 1 

Scalar transforms computation times 





TimeL= 128 


Time L = 256 


Time L = 512 


Time L= 1024 




(sec) 


(sec) 


(sec) 


(sec) 


Spharmonic- 


1.0e-01 


6.8e-01 


5.2e + 00 


4.2e+01 


Kit 


9.5e-02 


6.7e-01 


5.5e + 00 


4.6e+01 


HEALPix2.1 


4.8e-02 


3.0e-01 


2.0e + 00 


1.5e+01 


(i = 0) 


4.0e-02 


2.5e-01 


1.6e + 00 


l.le+01 


HEALPix2.1 


2.7e-01 


1.7e + 00 


l.le+01 


7.8e+01 


(1 = 1) 











NOTE. — Computation times for the SpharmonicKit and HEALPix2.1 
implementations of the scalar spherical harmonics transforms, measured on 
a 2.20 GHz Intel Pentium Xeon CPU with 2 Gb of RAM memory. Times 
associated with the direct transform are listed above corresponding times for 
the inverse transform. The HEALPix2. 1 direct transforms are reported both 
at (' = and i = 1 iteration. 



TABLE 2 

Scalar transforms error analysis 





Error L= 128 


Error L = 256 


Error L = 512 


Error L = 1024 


Spharmonic- 


3.3e-ll 


9.4e-ll 


2.6e-10 


7.5e-10 


Kit 


2.9e-09 


3.7e-09 


l.le-08 


1.8e-07 


HEALPix2.1 


7.0e-03 


4.8e-03 


3.5e-03 


2.4e-03 


(,• = 0) 


7.7e+00 


1.0e+01 


2.9e + 01 


2.3e + 01 


HEALPix2.1 


2.7e-04 


1.8e-04 


1.4e-04 


9.2e-05 


(1 = 1) 


2.7e-01 


3.6e-01 


l.le + 00 


7.7e-01 



NOTE. — Errors for the SpharmonicKit and HEALPix2.1 implementa- 
tions of the scalar spherical harmonics transforms, measured on a 2.20 GHz 
Intel Pentium Xeon CPU with 2 Gb of RAM memory. Relative root mean 
square errors after inverse and direct transforms are listed above the corre- 
sponding relative maximum errors on the coefficients. HEALPix2.1 results 
are reported both at i = and i = 1 iteration. 

up to L = 1024. But the relative root mean square errors are 
bounded by 7.0 x 10~ 3 , up to L = 1024, which ensures the 
global numerical precision of the implementation. This al- 
ready illustrates the numerical stability of the HEALPix im- 
plementation at the 0.7 % level, which however remains a 
much lower precision than with the SpharmonicKit imple- 
mentation. The use of the iterative scheme for the computa- 
tion of the direct transform allows to reach a better accuracy, 
at the price of a bigger time consumption. For the HEALPix 
implementation with one iteration (i =1), the relative maxi- 
mum errors on the coefficients and relative root mean square 
errors are then bounded by 1 . 1 and 2.7 x 10~ 4 respectively, up 
to L = 1024, therefore achieving a numerical stability at the 
0.027 % level. 

Let us emphasize that the conclusions of this comparison of 
memory requirements, computation times, and numerical sta- 
bility remain for the equi-angular and HEALPix implementa- 
tions of our standard correlation algorithm. Indeed, the over- 
all computation is strictly driven by the SpharmonicKit and 
HEALPix implementations of the scalar spherical harmonics 
transforms respectively. 

We finally briefly comment on the issue of changing the 
pixelization on which the signal is defined, from HEALPix 
to equi-angular, or conversely. A specific choice of pixeliza- 
tion might be preferred for some analyses. The above results 
confirm the very good precision with which both implementa- 
tions, with resolution N S ide = L/2 for HEALPix and on 2L x 2L 
equi-angular grids for SpharmonicKit, can calculate the direct 
and inverse transforms of a signal up to the band limit L. Con- 
sequently, a change of pixelization may be performed safely 
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by first computing the spherical harmonics coefficients of the 
signal by direct transform from one pixelization, and second 
operate an inverse transform on the alternative pixelization. 
The numerical precision of such an operation is obviously 
limited by the transform with the lowest precision, i.e. the 
HEALPix transform. 

Notice however that any spectral or directional correlation 
analysis only uses the spherical harmonics coefficients of the 
signal. In that perspective, the question of pixelization change 
makes no sense. The initial pixelization on which the signal is 
defined does not constrain such analyses as both the HEALPix 
and SpharmonicKit implementations allow a precise calcula- 
tion of the spherical harmonics coefficients of the signal. But 
the better numerical precision and lower computation times 
of the SpharmonicKit transform can really materialize if the 
data are originally sampled on equi-angular grids. Such a pos- 
sibility must be evaluated for each application independently. 
The question can notably be raised for CMB experiments such 
as WMAP and Pl anck, which currently use HEALPix grids 
JGorski et.al.200a. 

4.5. Existing alternative procedure and optimization with 
steerable filters 

For completeness of the discussion, we briefly expose 
an existing alternative algorithm for the directional corre- 
lation, based on the factorization of the rotati on operators 
R(p) on functions in L 2 (S 2 .dft) on the sphere as (Ris bo 19961 
iWandelt&Gorski 20011) 

r{ Vo A,x)=r(w-\ >-5>ft>)*(°>5>x+f)- (40) 

The directional correlation relation dl4> and the expression (|8) 
of the Wigner D-functions, matrix elements of the operators 
R(p), therefore give an alternative expression for the direc- 
tional correlation of arbitrary signals F and (non-steerable) 
filters on the sphere. We get indeed 

(R(p)V\F)= (^) mm ^ i(mVo+m ' eo+ " X \ (41) 

m,m' 

with the Fourier coefficients on the three-torus {R^\F) mm , n 
given by 

1>C 

(42) 

where C = max(|m|, \ m'\, \n\), and with the sym metry relation 
d' m ,J6) = d' mm ,(-9) (IVarshalovichetal. 19891) . For a band- 
limited signal F with band limit L on the sphere one has 
|m|,|m'|,|n| </ <L. In these terms, the directional correlation 
algorithm implemented through the factorization of rotations 
is structured as follows, with related asymptotic complexities: 

• Direct scalar spherical harmonics transforms, \l//„ and 
Fi m : 0(L?) on all pixelizations (0(L 2 logj L) on an equi- 
angular grid if the associated Legendre polynomials are 
pre-calculated). 

• Correlation {R^\F) mm , n in spherical harmonics space: 
C(L 4 ) through §2j. 

• Inverse Fourier transform (R^>\F) on the three-torus in 
(6>o, ifo): Oil? log 2 L) on an equi-angular grid by the use 
of the standard Cooley-Tukey fast Fourier transform al- 
gorithm, while 0(L 4 ) on HEALPix and other pixeliza- 
tions. 



The overall asymptotic complexity is therefore of order 0(L ) 
independently of the pixelization, due to the correlation in 
spherical harmonics space. 

Let us reconsider our first naive estimation of the compu- 
tation times for the directional correlation on maps of several 
megapixels on the sphere. For L ~ 10 3 , an 0(L 4 ) algorithm 
already greatly reduces computation times from years down 
to days on a single standard computer. However, if a large 
number of signals have to be considered, say thousands of 
simulations, this remains hardly affordable even through the 
use of multiple computers. 

Our algorithm achieves an 0(L 3 ) asymptotic complexity 
for the standard or directional correlation with steerable fil- 
ters, and even an OiL 2 logj L) asymptotic complexity on equi- 
angular grids if the associated Legendre polynomials are pre- 
calculated. The algorithm described here using the factoriza- 
tion of rotations can also achieve an 0(L 3 ) complexity for 
the standard or directional correlation, if steerable filters are 
considered. Indeed, for a steerable filter \& with an azimuthal 
band limit N, the index n decouples from asymptotic com- 
plexity counts: \n\ < N << L . But the corresponding struc- 
ture does not allow any reduction to 0(L 2 logjL), as it clearly 
appears from relation J42I . 

Notice that the first proposal for the d irectional 
wavelet analysis o f the CMB (McE wen et al. 2 005a: 
McE wen et al. 2 005b ) was based on the presently discussed 
algorithm using the factorization of rotations. Instead of 
considering steerable wavelet filters, a much less precise 
sampling in P points in the direction x was defined, with 
P << L, typically reducing the precision in the directional 
analysis from A\ = 2tt/L to Ax = 2ir/P. This also artifi- 
cially applies a band limit in the related Fourier index n: 
\n\ <P. Technically, the asymptotic complexity of directional 
correlation is thus reduced like if a steerable filter with an 
azimuthal band limit P was considered. However, this is 
by no means an answer to the problem of the directional 
correlation of band-limited signals on the sphere if one 
requires the same high precision (P ~ L) in the analysis of 
local directions as for the position of features on the sphere. 
In that perspecti ve, the use of stee rable filters is essential. 
As discussed in (Wiauxetal. 2005), the theoretical angular 
resolution power of a steerable filter is infinite in the sense 
that, around each point of the signal, it can resolve exactly the 
direction of local features with a pre-defined morphology. It 
therefore makes complete sense to define a precise sampling 
in the analysis of local directions (P ~ L). If needed, the 
explicit calculation of the directional correlation for all 
sampled values of \ therefore requires (see relation (12 1\ ) the 
add itional 0(P x L 2 ) ~ C(L 3 ) operation quoted in subsection 
14.31 Obviously, the resolution of the angular morphology 
itself of a local feature is limited by the non-zero angular 
width of the steerable filter, i.e. its angular band limit N in 
the index n associated with the azimuthal angle (p when the 
filter is located at the North pole. 

5. APPLICATION TO THE CMB SCALE-SPACE ANALYSIS 

5.1. Temperature and polarization 

The CMB radiation is completely described 
in terms of the observable temperature T and 
linear pol arization Q and U Stokes parameters 
ISeliak & Zaldarriaga 19971 IZaldarriaga & Seliak 1 997: 

IKosowskv 1996| [K amionkowski et al. 1997a; 

Kamionkowski et al. 1997bl iHu & White 19971 On the 
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one hand, the temperature is a scalar function on the sphere, 
i.e. invariant under local rotations in the plane tangent to 
the sphere at each point. It is also invariant under global 
inversion of the three-dimensional coordinates. On the other 
hand, the Q and U Stokes parameters are not invariant, but 
transform as the components of a transverse, symmetric, and 
traceless rank 2 tensor on the sphere, and respectively have 
even and odd parities under global inversion. Stated in other 
words, the combinations Q ± ill are spin ±2 functions on the 
sphere, and transform in one another under global inversion. 
In this context, the directional correlation of the temperature 
field with a scalar filter is invariant under the considered 
coordinates transformations. The directional correlations of 
the linear polarization parameters Q and U with a scalar filter 
are obviously not invariant, but transform in one another just 
as Q and U. A physical analysis will require the definition 
of invariant quantities. In that respect, the linear polarization 
of the CMB is equivalently defined by its electric E and 
magnetic B components, uniquely defined in terms of Q and 
U . These two polarization components are scalar functions 
on the sphere and also respectively have even and odd parities 
under global inversion. In this context, the statistics of a 
Gaussian and stationary CMB signal is completely described 
by the four invariant temperature (TT), polarization {EE and 
BE), and cross-correlation (TE) angular power spectra. In the 
same idea, the complete scale-space analysis of the CMB data 
through directional correlation on the sphere will therefore 
rely the directional correlation of the scalar temperature T, 
and the scalar polarization components E and B, with an 
arbitrary scalar filter. 

From the numerical point of view, E and B are related to the 
observables Q and U through derivative relations. It is con- 
sequently not possible to transform the experimental Q and 
U maps in scalar E and B maps from their basic definitions. 
However, the scalar spherical harmonics coefficients E/ m and 
B[ m are explicitly given as simple linear combinations of the 

spin ±2 spherical harmonics coefficients ±2(Q±iU) lm . From 
the recurrence relation ( I39> . these spin n = ±2 spherical har- 
monics coefficients may computed as linear combinations of 
\n\ + 1 = 3 scalar spherical harmonics coefficients. The direct 
scalar spherical harmonics transforms of E and B can there- 
fore be obtained from the experimental Q and U maps with 
the same asymptotic complexity as the direct scalar spherical 
harmonics transform of an experimental T map. From that 
stage, the remaining part of the detailed algorithmic structure 
proposed in § 14.31 can be applied identically to E[ m , B[ m , or 
7} m . In summary, our algorithm can be applied both for the 
directional correlation of temperature and polarization CMB 
maps with steerable filters. The overall asymptotic complex- 
ity for the analysis of polarization maps remains 0(L y ) on all 
pixelizations {0{L 2 \o^L) on an equi-angular grid if the as- 
sociated Legendre polynomials are pre-calculated). 

5.2. Numerical illustration 

As an illustration we apply our algorithm to the compu- 
tation of the wavelet coefficients of a temperature map with 
a second Gaussian derivative wavelet. A CMB temperature 
map is simulated both in the equi-angular and HEALPix pix- 
elizations, from the angular power spectrum which best fits 
the three-year WMAP data. A comparison of the implemen- 
tations of our algorithm on the two grids is proposed both in 
terms of computation times, and precision of analysis. 

First, we produce a simulated CMB temperature map both 



on equi-angular and HEALPix pixelizations in the following 
way. We define scalar spherical harmonics coefficients 7} m 
from a Gaussian distribution with zero mean and a variance 
given by the temperature angular po wer spectrum Cj T w hich 
best fits the three-year WMAP data ( Sperge leTal. 20061) . up 
to a band limit L = 1024. The spectrum Cj T was computed 
with CMBFAST 4 . Two real CMB temperature maps are then 
produced by inverse scalar spherical harmonics transforms 
of these coefficients, on a 2L x 2L equi-angular grid with 
the SpharmonicKit, and with a resolution N S ide = L/2 for the 
HEALPix pixelization. 

Second, we apply our algorithm to produce the directional 
correlation of each map with a second Gaussian derivative 
wavelet. The equi-angular or HEALPix implementations of 
our algorithm are used independently on the two correspond- 
ing maps. For coherence with the studied scalar spherical har- 
monics transforms codes, the complete equi-angular imple- 
mentation is coded in C, while the HEALPix implementation 
is coded in Fortran90. 

For a given scale adl*, the angular size of our dilated 
wavelets ^> a on the sphere is defined as twice the dispersion of 
the corresponding Gaussian. Fifteen scales are selected cor- 
responding to angular sizes of the second Gaussian deriva- 
tive lying between 20 arcminutes and 60 degrees: a\ — > 20', 
a 2 -> 40', a-, -> 1°, a 4 -► 5°, a 5 -> 10°, a 6 15°, a 7 -► 20°, 
a 8 25°, a 9 -► 30° am -► 35°, a u -► 40°, a n -> 45°, 
«i3 — > 50°, aw — > 55°, fli5 — > 60°. The analytical expression 
of the second Gaussian d erivative wavelet as a function of the 
scale a may be found in ( Wiaux et al. 2005). From a practi- 
cal point of view, the wavelet is first sampled on the sphere in 
each pixelization, and the corresponding spherical harmonics 

coefficients (^ a )im are calculated by direct spherical harmon- 
ics transform at each scale. 

The wavelet coefficients resulting from the directional cor- 
relation of the signal with the second Gaussian derivative 
at each scale a read as W^,{x,^o,a) = (/?(p) , 5 fl |r), for p = 
((po,9o,X), and ojq = ((po,0o). The coefficients resulting from 
the standard correlations with the corresponding basis fil- 
ters defined in i20\ read W^, (oJo, a ) = (R(wo)^ma\T), for 

1 < m < 3, and with * la = m 2a = tf"", and 

^3a = a * 9il<£auss \ The steerability relation i2\\ therefore ex- 
plicitly reads 

3 

Wl (x, wo, a) = k m (X) Wl m (fJo,a) , (43) 

m=l 

with the weights ki(x) = cos 2 x, feCx) = sm2 X' an d fe(x) = 
sin2% given in dl9l . At each of the fifteen scales considered, 
the coefficients w£ (u>o,a) gather all the information neces- 
sary to compute the directional correlation of the signal with 
the second Gaussian derivative. These coefficients are com- 
puted as functions on the sphere, up to the band limit L = 1024. 

On the one hand, the overall computation times at each 
scale for the coefficients (uo,<i) of the three basis filters 
(1 < m < 3) are of the order of 7.1 x 10 2 seconds in the equi- 
angular implementation on a 2.20 GHz Intel Pentium Xeon 
CPU with 2 Gb of RAM memory. The same calculation takes 
2.1 x 10 2 seconds in the HEALPix implementation wih zero 
iteration = 0) for the direct transforms. These values may 
directly be inferred from the computation times reported in 

4 See cmbfast.org for documentation (CMBFAST4.5.1 software). 
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Table ^ for the scalar spherical harmonics transforms. In- 
deed, the second Gaussian derivative contains the frequencies 
n = {0, ±2} and \n\ + 1 inverse scalar spherical harmonics are 
required for each spin- weighted transform in ( I29l l. Let us also 
recall that opposite spins are grouped by the algorithm before 
applying the scalar spherical harmonics transforms. For each 
of the three basis filters, four inverse scalar transforms (one 
at spin and three at spin 2) therefore add up to the direct 
scalar transform of the filter. Also adding the direct scalar 
transform of the signal, a total of four direct scalar transforms 
and twelve inverse scalar transforms are required, which leads 
to a very good estimation of the overall computation times re- 
ported above. 

On the other hand, the relative precision of the equi-angular 
and HEALPix implementations may be probed quantitatively 
through the calculation of the mean and variance of the 
wavelet coefficients at each scale a and each direction x, re- 
spectively 

[al]\ X ,a)=^-J^dn \W^(x,u J o,a)-^l(x,a)\ 2 .(44-) 

From relation J43> . those moments are explicitly given from 
the means and (co-)variances of the standard correlations at 
each scale a, respectively, as 

3 

M* (X, a) = ^2 k m (X) («) 

m=l 

3 

[4] 2 (X,a)= k m(x)k„Ax)vL>(a), (45) 

m,m' = l 

with 

Mm («) = T" / dQoWl (ujq , a) 
4?r J s i 

a Ln> ( fl )= ^ J d ^o (Wy m (w ,a)-/4„, («))* 

(V£ m ,(w ,a)-M* m >))- (46) 

All these quantities are defined from the directional corre- 
lation of a real signal with a real wavelet, and are therefore 
themselves real. Notice that quadrature rules need to be ap- 
plied in order to compute the latter integrals on both pixeliza- 
tions. On equi-angular grids, the sampling theorem discussed 
in |3 equivalently states that the integral of a scalar func- 
tion on the sphere with a band limit 2L may be computed 
exactly as a weighted sum of its sampled values on 2L x 2L 
equi-angular grids. The quadrature weights are the same as 
those appl ied for the computation of spherical harmonics co- 
efficients (Dris coll & Healv 19941) . On HEALPix grids, the 
integral of a function on the sphere is simply approximated 
by the sum of its sampled values, as all pixels have the same 
area. These rules are applied for the computation of the inte- 
grals /i m (a) an d ^mm' ( fl ) on tne sphere, at each of the fifteen 
scales considered. 

Figure Q] represents the six spectra <7 T mml (ci) computed in 
the equi-angular and HEALPix implementations, as functions 
of the scale. The top panel represents the three variances 
o\ lm (a), with 1 < m < 3. The bottom panel represents the 
three covariances o T mm ,(d) with 1 < m ^m! < 3. For any of 
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FIG. 1. — Comparison of the relative precision of the equi-angular and 
HEALPix implementations of the fast directional correlation algorithm. The 
six spectra associated with the wavelet coefficients from the directional cor- 
relation of a simulated CMB temperature map with the second Gaussian 
derivative wavelet are shown. The spectra, normalized as a' 2 a T ,(a) with 

r mm v * 

1 < m,m' < 3, are given in fiK 1 as functions of the wavelet scale a, con- 
verted in terms of angular size in degrees. The top panel represents the vari- 
ances cr^(a), o-JjCa), and cr^(a). The bottom panel represents the covari- 
ances aj 2 (a), and <rlL(a), For all six spectra and at any of the fifteen 
scales considered, the maximum relative difference between the two imple- 
mentations is bounded by 1.9 %. The equi-angular and HEALPix curves are 
therefore indistinguishable. 



those six spectra, at the first scale a\ — > 20', the maximum 
relative difference between the values arising from the equi- 
angular and HEALPix implementations is bounded by 1.9 
%. At that small scale, the assumption that the wavelet fil- 
ter is band-limited at L is not perfect anymore. The wavelet 
is not perfectly well sampled on the grids considered. Con- 
sequently, aliasing errors can occur in both the equi-angular 
and HEALPix implementations, which explain the discrep- 
ancy between the spectra. For the remaining fourteen scales, 
the maximum relative difference between the spectra from 
the equi-angular and HEALPix implementations drops below 
0.07 %. This remaining difference may be attributed to the 
calculation of the scalar spherical harmonics transforms in 
the HEALPix implementation, which is approximate even for 
band-limited signals and filters. This discrepancy is indeed 
coherent with the relative root mean square error on the spher- 
ical harmonics transform, estimated at 0.24 % in Table for 
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HEALPix with zero iteration (i = 0), and at L = 1024. We also 
know that this already very good precision of the HEALPix 
implementation can be enhanced using at least one iteration 
(i > 1), at the price of an increased computation time. This 
would reduce the relative differences with the theoretically 
exact equi-angular implementation. 

6. CONCLUSION 

We introduced a fast algorithm for the directional correla- 
tion of band-limited signals with band-limited steerable filters 
on the sphere. The a priori asymptotic complexity associ- 
ated with the directional correlation is 0(L 5 ), where 2L stands 
for the square-root of the number of sampling points on the 
sphere, also setting a band limit L for the signals and filters 
considered. The use of steerable filters allows to compute the 
directional correlation uniquely in terms of direct and inverse 
scalar spherical harmonics transforms. The overall asymp- 
totic complexity is correspondingly reduced to the asymptotic 
complexity of scalar spherical harmonics transforms, that is 
0(L?). We also emphasized that an 0(L 2 log\L) asymptotic 
complexity may be achieved on equi-angular pixelizations if 
associated Legendre polynomials are pre-calculated. We thor- 
oughly compared the implementations of the scalar spherical 
harmonics transforms on HEALPix and equi-angular grids, 
in terms of memory requirements, numerical stability, and 
computation times. First, the memory requirements are eas- 
ily accessible. Second, the SpharmonicKit implementation 
on equi-angular grids based on an exact algorithm proposed 
by Driscoll and Healy reaches the computer's numerical pre- 
cision. The exactness of the calculation relies on a sam- 
pling theorem for scalar functions on equi-angular grids on 
the sphere. The already good precision of the HEALPix im- 
plementation of the scalar spherical harmonics transform may 
be enhanced through an iterative process, but never reaches 
the level achieved on equi-angular grids. Third, the compu- 
tation times for the 0{L?) scalar transform of maps of sev- 
eral megapixels on the sphere (L ~ 10 3 ) are reduced from 
years to tens of seconds on a single standard computer in 



both the equi-angular and HEALPix implementations. The 
equi-angular implementation is slightly less rapid than the 
HEALPix implementation with zero iteration (i = 0), but be- 
comes more rapid as soon as at least one iteration (i > 1) is 
considered in the HEALPix scheme. The computation of the 
directional correlation of multiple signals or simulations with 
steerable filters is thereby rendered easily affordable at high 
resolutions. 

In the perspective of the scale-space analysis of the CMB 
temperature (T) and polarization (E and E) anisotropies, the 
wavelet processing of the CMB data is therefore easily afford- 
able at the high resolutions of the WMAP and Planck exper- 
iments. The identification of possible local non-Gaussianity 
or statistical anisotropy signatures, or foreground emissions, 
is accessible with the same high precision in both position 
and local direction on the sphere. The low computation times 
and very good precision of the equi-angular and HEALPix 
implementations of our algorithm was clearly illustrated in 
that context through a wavelet decomposition of a simulated 
three-year WMAP temperature map of several megapixels. 

The generic algorithm developed for the fast directional 
correlation on the sphere may also find other applications, no- 
tably in the analysis of asymmetric beam effects on the CMB 
temperature and polarization data, or well beyond cosmology. 
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